arXiv:1508.01504vl [cs.DS] 6 Aug 2015 


Resource Oblivious Sorting on Multicores * 


Richard Cole 1 Vijaya Ramachandran * 

August 7, 2015 


Abstract 

We present a deterministic sorting algorithm, SPMS (Sample, Partition, and Merge Sort), 
that interleaves the partitioning of a sample sort with merging. Sequentially, it sorts n elements 
in O(nlogn) time cache-obliviously with an optimal number of cache misses. The parallel 
complexity (or critical path length) of the algorithm is 0(logn • log log n), which improves on 
previous bounds for optimal cache oblivious sorting. The algorithm also has low false sharing 
costs. When scheduled by a work-stealing scheduler in a multicore computing environment with 
a global shared memory and p cores, each having a cache of size M organized in blocks of size 
B, the costs of the additional cache misses and false sharing misses due to this parallel execution 
are bounded by the cost of 0(S ■ M/B) and 0(S ■ B) cache misses respectively, where S is the 
number of steals performed during the execution. Finally, SPMS is resource oblivious in that 
the dependence on machine parameters appear only in the analysis of its performance, and not 
within the algorithm itself. 


1 Introduction 

We present a parallel sorting algorithm, which we call Sample, Partition, and Merge Sort (SPMS). 
It has a critical path length of O(lognloglogn) and performs optimal 0{n log n) operations with 
optimal sequential cache misses. Further, it has low caching and false sharing overhead when 
scheduled by a work stealing scheduler on a multicore with private caches, and is resource oblivious, 
i.e., it uses no machine parameters in its specification. 

At the heart of the sorting algorithm is a recursive multi-way merging procedure. A notable 
and novel aspect of this procedure is that it creates its recursive subproblems using a sample sort 
methodology. We view SPMS as interleaving a merge sort with a sample sort in a natural way. 

Previous Work. Sorting is a fundamental algorithmic problem, and has been studied extensively. 
For our purposes, the most relevant results are sequential cache-oblivious sorting, for which provably 
optimal algorithms are given in P2!> optimal sorting algorithms addressing pure parallelism mm, 
and recent work on multicore sorting mmmm- When considering the parallelism in an algorithm, 
we will interchangeably use the terms critical pathlength, span, and parallel time to denote the 

*A preliminary version of this paper appeared in Proc. International Colloquium of Automata, Languages and 
Programming (ICALP), Track A, Springer LNCS Volume 6198, pp. 226-237, 2010. 

'Computer Science Dept., Courant Institute of Mathematical Sciences, NYU, New York, NY 10012. Email: 
cole@cs.nyu.edu. This work was supported in part by NSF Grants CCF-0830516 and CCF-1217989. 

^Dept. of Computer Sciences, University of Texas, Austin, TX 78712. Email: vlr@cs.utexas.edu. This work 
was supported in part by NSF Grants CCF-0830737 and CCF-1320675. 


1 



number of parallel steps in the algorithm when it is executed with a sufficiently large number of 
processors to exploit all of the parallelism present in it. 

The existing multicore algorithms take two main approaches. The first is merge sort H 01, 
either simple or the pipelined method from HDJ. The second is deterministic sampling J22J [6]: 
this approach splits the input into subsets, sorts the subsets, samples the sorted subsets, sorts the 
sample, partitions about a subsample, and recursively sorts the resulting sets. Our algorithm can 
be viewed as applying this approach to the problem of merging a suitable number of sorted sets, 
which eliminates the need for the first two steps, resulting in a smaller critical pathlength. 

More specifically, the algorithm in [3] is a simple multicore mergesort; it has polylog parallel 
time, and good, though not optimal cache efficiency; it is cache-oblivious for private caches (the 
model we consider in this paper). The algorithm in [3] achieves the optimal caching bound on 
an input of length n, with O(logn) parallel time (modulo dependence on cache parameters), but 
it is both cache-aware and core-aware; this algorithm is based on m- The algorithm in [22] is 
designed for a BSP-style version of a cache aware, multi-level multicore. It uses a different collection 
of parameters, and so it is difficult to compare with it directly. The algorithm in [6j recursively 
sorts the samples, achieving 0(log 2 n) critical path length deterministically, and 0(log 1,5 n) in a 
randomized version. 

Our Results. Our main results are stated in the following two theorems. Theorem 11.11 which 
applies to any work-stealing scheduler, is established in Section [6] and Theorem ll.21 which applies 
to the randomized work-stealing scheduler (RWS), is established in Section [7] These theorems use 
the following notation: M is the cache size, B is the block size, a steal is the process of transferrring 
a parallel task from the processor that generated it to another processor that will start its execution, 
p is the number of processors, b is the cost of a cache miss, and s > b is the cost of a steal. 

Theorem 1.1. There is an implementation of the SPMS algorithm that performs 0(n log n) oper¬ 
ations, runs in parallel time O(lognloglogn), and, under work-stealing, incurs 0(j| iogM § -s) 
cache misses, and incurs a total cost due to false sharing bounded by the cost of 0(S■ B) = 0{^ ■ S) 
cache misses, where S is the number of steals. These results assume a tall cache, i.e., M = Q(B 2 ). 

Theorem 1.2. W.h.p. in n, the execution time of SPMS with the RWS scheduler, when including 

£ log B \\ 
b ' B }) 

ifs = 0{b •§) and M = n(B 2 ). 

Roadmap. In Section [2] we describe the basic SPMS sorting algorithm, and show that it has a 
fairly simple implementation with optimal work (i.e., sequential running time) and O(logradoglogn) 
critical pathlength, and also a variant with optimal work and optimal sequential cache complexity. 
These two sets of bounds are obtained with different implementations. In Section [4] we present an 
implementation that achieves all three bounds simultaneously. This section follows Section [3] which 
describes the multicore caching model, together with some basic methodology used in prior work for 
developing efficient multicore algorithms HU3UEU6UI3], and some basic framework from [13] relevant 
to the multicore caching model. Then, in Sections [5] and [6] we address the overhead of false sharing , 
a phenomenon that is inherent in a resource oblivious multicore setting: Section [5] bounds the false 
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sharing costs in SPMS when there is no deallocation of space (as in, e.g., a functional programming 
setting), and Section [6] bounds the cost of additional fs misses that occur on the execution stack 
when space is re-used there; the background for this setting is given in Section 13.31 The results 
in Sections [MS] build on our earlier work on false sharing in [13] for hierarchical balanced parallel 
(HBP) computations; however, since SPMS is not exactly an HBP computation due to variations 
in sub-problem sizes, we present a self-contained analysis of the false sharing results for SPMS. 
In Section [7] we analyze the performance of our final SPMS implementation (given in Section [5]) 
as a function of the number of parallel tasks scheduled by the scheduler. We also derive specific 
results for the Randomized Work Stealing (RWS) scheduler, using some already known performance 
results. We conclude in Section [8] with a summary, and an avenue for further research. 

2 SPMS, A New Deterministic Sample, Partition, and Merge Sort 

There are two main approaches used in parallel algorithms for sorting (ignoring methods used in 
sorting networks): Sorting by merging, based on merge-sort, and sorting by sampling and partition¬ 
ing, called Sample Sort or Distribution Sort. Both methods can also achieve good cache-efficiency. 

In parallel merge sort, the input is grouped into two or more sets, which are sorted recursively, 
and the sorted sequences are then merged by a parallel merging algorithm. If the merge is binary 
then we obtain a parallelization of merge-sort, which can achieve optimal O(logn) parallel time 
[10], though it is not known how to achieve optimal cache-oblivious cache-efficiency. Multi-way 
merge can achieve optimal cache-oblivious cache-efficiency, but it is not known whether very high 
levels of parallelism can be achieved when cache-oblivious cache-efficiency is also desired. 

In sample sort, a small sample of the input is taken and then sorted, typically by a simple 
parallel algorithm, and not recursively since the sample size is small relative to the input size. The 
sorted sample elements are used to partition the input elements into a linearly ordered sequence 
of subsets of the input elements, and finally each of these subsets is sorted recursively. Sample (or 
distribution) sort can achieve optimal cache-oblivious cache-efficiency. 

In sample sort, all of the recursively sorted subsets have approximately the same size, and this 
is achieved by having the ranks of the samples spaced fairly evenly among the input elements. This 
is readily ensured by choosing a random sample of the desired size. 

In deterministic sample sort, the sample is usually chosen by partitioning the input into a 
collection of sets, sorting each set (recursively), then selecting every £;-th element (for a suitable 
k) as a sample element, and then sorting the resulting samples from the different sets by a simple 
algorithm. Often, ‘oversampling’ is employed, where every Z-th element in the sorted sample is used 
in the actual sample. Oversampling with suitable parameters ensures that the partitioned subsets 
will have approximately the same size. 

As noted earlier, after the sorted sample is obtained, the input is partitioned into a linearly 
ordered sequence of subsets, each of which is sorted recursively. However, in deterministic sample 
sort, this is somewhat wasteful, since each set (from which the sample elements were chosen) was 
sorted in order to obtain evenly spaced sample elements, yet the sorted information in these sets is 
not used once the samples are selected. 

Our algorithm, Sample, Partition and Merge Sort (SPMS), uses a variant of the sample and 
partition method of deterministic sample sort, but it maintains the sorted information from each 
original sorted subsequence within each linearly ordered subset. Hence, the recursive problem to 
be solved is not a recursive sort, but a recursive multi-way merge. Accordingly, we generalize the 
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problem being solved as follows. Its input A comprises up to r sorted lists of total length n < 3 • r c 
where c > 6 is a constant, and the task is to merge these lists. We let SPMS (A; n,r) denote this 
problem instance. A sort is executed by setting r = n, i.e. the input comprises n lists each of length 
1. 

The SPMS algorithm performs two successive collections of recursive y/r- way merges, each 
merge being of size at most 3 • W 2 . To enable this, suitable samples of the input lists will be sorted 
by a logarithmic time procedure, which then allows the original problem to be partitioned into 
smaller subproblems that are merged recursively. Here is the high-level algorithm. 

SPMS(i;n,r) 

Input. A is a collection of r sorted lists L\, L 2 , ■ ■ ■ , L r , of total length n < 3 • r c . 

Output. The elements of A rearranged in sorted order. 

Step 0. if n < 24 then sort with any sequential 0{n log n) time algorithm else perform Steps 1-3. 

Step 1. If n > 3W 2 then Deterministic Sample and Partition: 

la. Sample. Form a sample S of every r( c / 2 ) -1 -th element in each of the sorted lists in A, and 
sort S. 

lb. Partition. Form a pivot set P containing every 2r-th element in the sorted S (and for 
convenience add a zero’th pivot with value —00 and a new largest pivot with value 00 ). Use P to 
partition the original problem A into k = 0(n/r 2 ) disjoint merging subproblems, A 2 , ■ ■ ■ , A/., 
each comprising at most r sorted lists, with the items in A* preceding those in Aj+i, for 1 < i < k. 
Thus, the items in Aj are those elements in A with values at least as large as the (i — l)st pivot, 
and less than the value of the i-th pivot. We show below that each Aj contains at most 3c • W 2 
elements. 

Step 2. First Recursive Multi-way Merge: 

For each subproblem Aj, group its lists into disjoint subsets of at most y/r lists, and merge the lists 
in each group. As Aj contains at most 3 • r^ items, this bound applies to each of the individual 
groups too. Thus the y/r- way merge in each group can be performed recursively. The output, for 
each subproblem Aj, is a collection of at most y/r sorted lists of total length at most 3 • n. Thus 
Step 2 is the following. 

2. for i = 1 to k do 

2a. Form disjoint A t j , with j < y/r, where each Ajj contains at most y/r of the lists in Aj; 
let |Ajj| = rijj. for each j. 

2b. for each subproblem Ajj do SPMS ( Ajj ; njj, y/r) 

Step 3. Second Recursive Multi-way Merge. 

for each subproblem Aj do recursively merge the y/r sorted lists computed in Step 2. 

Let | Aj | = n t for each i. Step 3 is the following. 

3. for i = 1 to k do SPMS(Aj; m, y/r). 

Comment. The subproblems in Step 2 may have quite different sizes, for these depend on the exact 
positions of the sampled items, which the algorithm cannot completely control. This, incidentally, 
is the reason for focussing on the more general merging problem with possibly unequal sized sorted 
lists. 

Lemma 2.1. Every subproblem Aj created in Step lb has size between r c / 2 + 1 and 3 • W 2 — 1, 
except possibly the last subproblem, which may be smaller. 
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Proof. Let S be the subset of S comprising every 2r-th item in S used in Step lb. Let e and e! 
be two successive elements in sorted order in S. We will argue that their ranks in the full set of 
n items differ by between ?’ c / 2 + 1 and 3 • W 2 — 1. Consider the j-th input list Lj and let SLj be 
the sorted list of samples in S from Lj. Suppose there are kj items from SLj in the range [e, e'\. 
Then there are between (kj — 1) • r c / 2_1 + 1 and (kj + 1) • W 2-1 — 1 items from Lj in this range. 
Summing over all j, as the sum of the kj ’s totals 2r, we obtain that there are more than r c / 2 + r 
items in this range and at most 3 • W 2 — r items in this range. The upper bound can be tightened 
by 2(r — 1) on noting that the upper bound is not met for the lists containing e and e'. ■ 

In the rest of the paper we give efficient implementations of SPMS for cache-oblivious cache 
efficiency, for work-efficient parallelism, for simultaneously achieving both of the above, and for 
reducing the communication cost of false sharing. In all of these implementations, Steps 2 and 
3 remain unchanged, but we give several different implementations of Step 1 with the goal of 
achieving efficiency with the most direct method for the measure being considered. This includes 
several different methods for sorting the sample S efficiently in Step la. 

2.1 A Cache Oblivious Sequential Implementation 

We consider a standard caching environment as in m We assume that data is organized in blocks 
of size B, and an ideal cache of size M; here, an ideal cache is one that has an optimal cache 
replacement policy. We assume the standard tall cache assumption, namely that M = Ll(B 2 ), so 
the cache can hold Ll(B) blocks. Since, with a tall cache, n/B dominates yfn when n > M, and 
when n < M, a cache bound of n/B implies that each data item is read into cache only once, we 
have the following fact. 

Fact 2.2. With an ideal tall cache, a cache miss bound of n/B + y/n can be bounded as 0(n/B). 

As usual, we will refer to a cache miss bound of 0(n/B ) for accessing n elements as the scan 
bound and a cache miss bound of 0(^\og M n) as the sort bound; each of these two bounds is 
optimal for the corresponding problems of scanning and sorting. 

Data Layout. Data layout is important for cache-efficiency. We will assume that the r input lists 
are stored in an array, starting with the sorted items in L\ , followed by L2, and so on. Similarly, 
we will store the sample in an array S’fL.s], which will be a subsequence of the input array. 

We now describe an efficient cache-oblivious implementation of Steps la and lb. For Step la, 
we will compute the rank of each element in S, which is essentially equivalent to sorting S. We will 
not perform the final step of storing S in sorted order, so that we can avoid the cache misses that 
would result with this step. 

Step 1 of Cache-oblivious SPMS 

Step la. Form the sample S and rank each element of S in S. 

(i) For each of the r sorted lists Lj, form the sublist SLj comprising every r( c / 2 ) -1 th item in 
Lj. The set S comprises the union of the items in the SLj lists. 

(ii) For each item in S, determine its rank in each SLj, storing the ranks for each item con¬ 
tiguously as follows. 

Let |S| = s. Use an array of size r • s, where the first r locations will store the rank of £[1] 
in each SLj, 1 < j < r, and each subsequent group of r locations will store the r ranks of each 
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subsequent element in S. Scan S r times, once in turn for each of the lists SL \, SL 2 ,..., SL r ; the 
j-th scan, which scans S once and SLj r times, records the rank in SLj of each item e in S. 

(Hi) For each e € S, sum its r ranks computed in Step la (ii) to obtain its rank in set S. 

Step lb. Form the sorted pivot set P and partition the lists with respect to elements in P. 

( i ) Store the pivot set P in sorted order by scanning the ranks of the items in S, and writing 
those items with rank an integer multiple of 2r to their sorted location in an array of size s/(2r). 
In general, contiguous inputs can be written to non-contiguous output locations (and vice versa). 

(ii) Extract each subproblem A, in turn to complete the partitioning. This will incur some 
caching penalty over the scan bound because for each subproblem fragments of r input lists are 
read and written. 

The following observation will be used in deriving our cache miss bound for the above steps. 
Observation 2.3. If r > B and s = n/A c / 2 ) -1 ; then r ■ s = 0(n/B) if c > 6. 

Proof, r ■ s = r ■ n/r ( c / 2 ) -1 = 0(n/(B ■ A c / 2 ^~ 3 )) = 0(n/B) when c > 6. ■ 

Lemma 2.4. SPMS can be implemented to run in optimal sequential cache-oblivious cache com¬ 
plexity 0((n/B) ■ log M n) if c > 6 and M > B 2 . 

Proof. Recall that s = 151. 

Step la(i) is simply a scan of the full set of items and so incurs 0(n/B) cache misses. Step 
la(ra) is a segmented sum computation on an array of size sr and so can be performed in O(srfB) 
cache misses. If c > 4, s ■ r = n/A c t 2 > 2 , hence s ■ r / B = 0(n/B) when c > 4. 

Consider the reads in Step la (ii). This step performs r scans over the S array, and over each 
of the SLj for the reads. This has caching cost s ■ r/B + Ylj r ■ \SLj\/B. As already noted, 
s ■ r/B = 0(n/B) when c > 4. To bound the caching cost for accessing the SLj . note that the 
total cost for the r accesses to those SLj that have at least B items is 0(r ■ s/B) since s is the 
sum of the sizes of all of the SLj, and \SLj\/B > 1 for each of these lists. For any SLj with fewer 
than B items, we observe that only one cache miss is needed to bring this list into cache, and in an 
ideal cache, this list will stay in cache during the entire scan of S that determines the ranks of the 
elements of S in this SLj, as long as the number of such lists is less than M/B. Hence the cost for 
accessing this class of lists is just r cache misses, which is 0(\/n) when c > 4 since n = Ll(r c ^ 2 ). If 
r > M/B then the tall cache assumption implies r > B, and by Observation 12.31 r • s = 0(n/B) if 
c > 6. 

Now we bound the writes in Step la (ii). There are r • s writes in this step, and each sequence 
of s writes is written to positions r apart in the output array. If r > B, r ■ s = 0(n/B) if c > 6 by 
Observation 12.31 If r < B, then this step writes B/r successive items with one cache miss, hence 
the cache miss cost is bounded by 0(r ■ s ■ (r/B)) = 0(n/(r c / 2 ~ 3 B)), and this is 0(n/B) when 
c > 6. 

In Step lb(i) the reads have cost 0(n/B) since the array being scanned has size 0(n/rl c / 2 ) -1 ). 
Every 2r-th of the scanned elements is written, and as each write could incur a cache miss, there 
are 0(n/r c / 2 ) cache misses for the writes. This is 0(y/n) since n < r c . 

Finally, for Step lb(iz), as the subproblems are written in sequence, we can bound the cache 
misses as follows. If r < B, it incurs only \Ai\/B = 0(n/B) cache misses since M > B 2 and so 
the (ideal) cache can hold the most recently accessed block from each of the r < B lists. If r > B, 
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it incurs JT[|Aj|/.B + O(r)] = n/B + JT 0(r) cache misses. But if r > B, JT 0(r) = 0(^2 ' r ) = 
0(n/r c / 2 ~ 1 ) = 0(n/r) = 0(n/B), if c > 4. 

Thus for Step 1, we have a cache miss bound of 0(n/B + y/n), which is 0(n/B) by Fact 12.21 
This yields the following recurrence for Q(r,n), the number of cache misses while executing an 
instance of SPMS on an input of size n < r c (recall that SPMS calls the multi-way merge with 
n lists, each containing one element, and with n = r initially). In this recurrence is the size 
of the multi-way merge subproblem A t] in Step 2 , and n ? ; is the size of the z-th multi-way merge 
subproblem in Step 3. 


Q(n,r ) = 0(n/B) + ^Q(nij,y/r) + ^Q{ni,y/r). 

i,j i 

The base case, n < M, has Q(n,r) = 0(n/B). The solution to this recurrence is the desired 
optimal caching bound Q(n,r) = O (jj \og m ^J • ® 

Comment. The bound in Lemma [2.41 can be generalized to handle a smaller tall cache assumption, 
namely an assumption that M > B 1+€ for any given constant e > 0, as follows. We will now need 
that c > 2 + 2 • —. The analysis of Step 1 now splits into cases according as r is smaller or larger 
than B e . Observation 12.31 which is used in the analysis of Step la(ii), now holds if r > B e ] also 

in the analysis of Step lb(ii), we obtain that when r > B e , ^Tr = O ^^72 • ^ = 0(n/r 3_1 ) = 
0(n/r 1//e ) = 0(n/B). 

2.2 A Parallel Implementation 

Here we follow the framework of the cache-oblivious algorithm in the previous section. The imple¬ 
mentation is changed in a few places as follows. 

In Step la (ii) we perform r binary searches for each item e in SL to determine its rank in each 
of the lists SLj. To carry out Step lb(u), for each item in the input, we will compute its destination 
location in the partitioning; then it takes a further 0 ( 1 ) parallel time to perform the partitioning. 

To this end, we find the rank of each of the n/(2 • W 2 ) pivot items in each list Lj (by means of 
n /(2 • r C//2_1 ) binary searches) and we also find for each pivot item e, the destination rank of its 
successor in each list Lj. These can be computed using a prefix sums computation following the 
binary searches. Then a parallel scan across the lists Lj yields the destination locations. The rest 
of the computation in Steps la and lb is readily performed in logarithmic time and linear work. 
This leads to the following lemma. 

Lemma 2.5. There is an implementation of SPMS that uses 0[n log n) operations and 0(log n log log n) 
parallel steps on an input of size n, if c > 6. 

Proof. Clearly, the work, i.e., the cost or operation count, and parallel time for Step 1, apart from 
the binary searches, are O(n) and O(logn) respectively. Each binary search takes O(logr) time 
and there are 0(r ■ n/r c / 2-1 ) binary searches, giving an overall work of 0 (n/r c / 2 ” 2 ) for the binary 
searches; this is an 0(n ) cost if c > 4. 

Let W (r, n) be the operation count for a collection of multi-way merging problems of total size 
n, where each comprises the merge of at most r lists of combined size at most r c . Then we have: 
W(r,n ) < 0(n) + 2IF(r 1 / 2 ,n) = 0(n log r) = 0{n log n). 
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For the parallel time, we have seen that Step 1 runs in O(logn) time. Hence, the parallel run 
time T(r,n) is given by: T(r,n ) < O(logn) + 2T(y/r, W 2 ) = O(lognloglogr) = O(lognloglogn). 


Lemmas 12.41 and [231 lead to the following corollary. 

Corollary 2.6. 

i. If M = It(B 1+e ), for any given constant e > 0, SPMS can be implemented to sort n elements 
cache-obliviously in O(nlogn) time with 0((n/B) log M n) cache misses. 

ii. SPMS can be implemented to sort n elements in 0(log n • log log n) parallel steps and 0{n log n) 
work. 

In Section 0] we will describe an implementation of Step 1 that will give rise to a parallel, 
cache-oblivious version of SPMS that simultaneously attains the above work, parallel time, and 
cache-oblivious cache-efficiency (when considering only the standard caching costs considered in 
sequential analysis). Then, in the following sections, we describe and bound the effects of ‘false- 
sharing’. 

3 Multicore Computation Model 

We consider a multicore environment consisting of p cores (or processors), each with a private cache 
of size M. The p processors communicate through an arbitrarily large shared memory. Data is 
organized in blocks (or ‘cache lines’) of size B. 

If we ignore communication costs a multicore is simply an asynchronous PRAM [19] . In contrast 
to a completely asynchronous environment, in this paper we use a computation model in which 
synchronization is performed using binary forks and joins; we describe this model in Section 13.11 
Parallel tasks are assigned to idle processors by a scheduler, and we discuss the type of schedulers 
we consider in Section m The communication costs in the multicore arise from the delay when 
transferring a requested piece of data from shared memory to the cache of the processor requesting 
the data. We discuss the caching costs in Section l3Tl We discuss the notion of resource obliviousness 
in Section 13.51 Finally, in Section 13.61 we discuss related models and schedulers as well as the 
connections between this current paper and its earlier proceedings version in HU. 

3.1 Computation Dag Model 

We will express parallelism through paired fork and join operations. A fork spawns two tasks that 
can execute in parallel. Its corresponding join is a synchronization point: both of the spawned 
tasks must complete before the computation can proceed beyond this join. 

The building block for our algorithms, which we call a fork-join computation, comprises a 
height O(logn) binary tree of fork nodes with possibly 0(1) additional computation at each node, 
with n leaf nodes, each performing a sequential computation of length 0(1) or O(logn), followed 
by a complementary tree of join nodes on the same leaves, again with possibly 0(1) additional 
computation at each node. In general, one might have other lengths of sequential computation, 
but they do not occur in the SPMS algorithm. 

A simple example of a fork-join computation is the computation of the sum of the entries in an 
array, by means of a balanced tree, with one leaf node for each array entry, and each internal node 


returning to its parent the sum of the values computed by its children. Good cache performance 
occurs if, as is natural, the left-to-right ordering of the leaves matches the order of the entries in 
the array. Prefix sums can be computed by sequencing two fork-join computations. 

The overall SPMS algorithm is created by means of sequencing and recursion, with the fork-join 
computations forming the non-recursive portions of the computation. Such a computation contains 
only nested fork-join computations, and it generates a series parallel graph. The computation dag 
for a computation on a given input is the acyclic graph that results when we have a vertex (or node) 
for each unit (or constant) time computation, and a directed edge from a vertex u to a vertex v if 
vertex v can begin its computation only after u completes its computation. For an overview of this 
model, see Chapter 27 in fl5] , 

3.2 Scheduling Parallel Tasks 

During the execution of a computation dag on a given input, a parallel task is created each time a 
fork step / is executed. At this point the main computation proceeds with the left child of / (as in 
a standard sequential dfs computation), while the task r at the right child r of the fork node is a 
task available to be scheduled in parallel with the main computation. This parallel task r consists 
of all of the computation starting at r and ending at the step before the join corresponding to the 
fork step /. 

Consider the execution on p processors of a computation dag D on a given input. In general, 
the parallel tasks generated during this execution can be scheduled across the p processors in many 
different ways, depending on the policy of the scheduler used to schedule parallel tasks. For example, 
the task r mentioned in the above paragraph could be moved by the scheduler to a processor Q 
different from P , the processor executing the main computation, in which case t will execute on 
Q in parallel with the execution of the main computation on P. The scheduler could also choose 
not to schedule r on another processor (perhaps because all other processors are already executing 
other computations), and in this case, r will be executed by P according to the sequential execution 
order. In the case when the parallel task r is moved from P to another processor Q, we will say that 
processor Q steals r from P; this terminology is taken from the class of work-stealing schedulers, 
which is the class of schedulers we consider in this paper. 

Work-stealing Schedulers. Under work-stealing, each processor maintains a task queue on 
which it enqueues the parallel tasks it generates. When a processor is idle it attempts to obtain 
a parallel task from the head of the task queue of another processor that has generated parallel 
tasks. The exact method of identifying the processor from which to obtain an available parallel task 
determines the type of work-stealing scheduler being used. The most popular type is randomized 
work-stealing (RWS, see e.g., [7j), where a processor picks a random processor and steals the task 
at the head of its task queue, if there is one. Otherwise, it continues to pick random processors 
and tries to find an available parallel task until it succeeds, or the computation completes. 

3.3 Execution Stacks 

Let us consider how one stores variables that are generated during the execution of the algorithm. 
It is natural for the original task and each stolen task to each have an execution stack on which 
they store the variables declared by their residual task. E r will denote the execution stack for a 
task r. As is standard, each procedure and each fork node stores the variables it declares in a 
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segment on the current top of the execution stack for the task to which it belongs, following the 
usual mode for a procedural language. As implied by its name, an execution stack is accessed by 
a processor in stack order. 

Execution Stack and Task Queue. The parallel tasks for a processor P are enqueued on its 
task queue in the order in which their associated segments are created on P’s execution stack. 
The task queue is a double-end queue, and P will access it in stack order similar to its accesses to 
its execution stack. Thus, P will remove an enqueued task a from its task queue when it begins 
computing on cr’s segment. 

Under work-stealing, the task that is stolen (i.e., transferred to another processor) is always the 
one that is at the top of the task queue in the processor from which it is stolen, or equivalently, 
tasks are stolen from a task queue in queue order. 

Usurpations of an Execution Stack. Suppose that a subtask t' is stolen from a task r. 
Consider the join node v immediately following the node at which the computation of t' terminates. 
Let C be the processor executing t — t' when it reaches node v and let C be the processor executing 
t' at this point. To avoid unnecessary waiting, whichever processor (of C and C') reaches v second 
is the one that continues executing the remainder of r. If this processor is C' , we say that C' has 
usurped the computation of r. We will examine the additional false sharing costs incurred due to 
usurpations in Section [6l There are additional regular caching costs as well, but these are handled 
by the analysis described below in Section T3. 4.21 

3.4 Caching Cost 

As in a sequential computation, a processor sustains a cache miss when accessing a data item if 
that item is not in its cache. When we bound the cache miss cost of a parallel execution, we will 
assume an optimal cache replacement policy (as in the sequential case). The caching cost for a 
parallel execution is the sum of the caching costs at each of the p parallel cores. There is a second 
source of cache misses that is present only in the parallel setting: false sharing. We will discuss 
and analyze this cost in Sections [SHSJ Here we confine our attention to regular cache misses. 

3.4.1 Execution Stack and Cache Misses due to Misalignment 

A parallel execution can incur additional cache misses due to block misalignment. Block misalign¬ 
ment can arise whenever there is a new execution stack, for the block boundaries on the new stack 
may not match those in the sequential execution. In fact, block misalignment, which was overlooked 
previously, arises even in analyses in the style of [Tj. But, as we will see, in Corollary 13.21 the cost 
due to block misalignment is at most the cost increase that occurs if the cache sizes are reduced by 
a constant factor, and for regular algorithms m, which cover most standard algorithms including 
SPMS, this increases the cache miss costs by only a constant factor. 

Until now our analysis had assumed that the parallel execution created the same blocks as 
the sequential execution and argued that the resulting cache miss bound was of the form 0(Q + 
additional terms), where Q was the cache miss cost of a sequential execution. But, in fact, this 
assumption need not hold. 

Let us consider what happens when a subtask a is stolen from task r. The variables that a 
declares are placed on its execution stack, presumably starting at the boundary of a new block. 
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However, if there had been no steal, these variables would have been placed on the continuation of 
r’s execution stack, presumably, in general, not starting at a block boundary. This misalignment 
may cause additional cache misses compared to the sequential execution. 

More specifically, any block /3 on the execution stack in the sequential algorithm may be ap¬ 
portioned to multiple execution stacks in a parallel execution, with one of these stacks potentially 
using portions of two of its blocks to hold its portion of /3. 

To bound the cost of cache misses, it will suffice to bound, for the processor P a executing a , 
for each block /3 in the sequential execution storing data that P a needs to access, the number of 
blocks storing portions of (3 that P a accesses in the parallel execution. 

Lemma 3.1. Let P a be a processor executing a stolen subtask a. Suppose P a accesses data which 
is stored in block /3 in the sequential execution. Then P a needs to access at most 4 blocks in the 
current parallel execution to access this data. 

Proof. To analyze this, let us consider which variables P a can access. P a can access the 0(1) 
variables declared by the node from which a was forked. It can also access variables declared by 
its “calling” SPMS task (the SPMS task that owns the steal that created a). Finally, it can access 
the variables stored on cr’s execution stack. These variables could be on as many as three different 
execution stacks, one for the calling task, one for the parent fork node, and one for a. 

Consider a block /3 in the sequential execution, some of whose contents P a seeks to access. If /3’s 
contents are stored in a single execution stack in the current parallel execution, then /3’s contents 
occupy at most two blocks in the current parallel execution. On the other hand, if in the current 
parallel execution, /3’s contents are stored over two of the execution stacks that P a accesses, then 
/3’s contents occupy at most three blocks that P a accesses (this can happen for at most two distinct 
blocks /3 accessed by P a ). Finally, if /3’s contents are stored over three execution stacks that P a 
accesses, then /3’s contents occupy at most four blocks that P a accesses (this can happen only once 
for data being accessed by P a , and only if the “over two execution stacks” scenario does not occur). 

It follows that for each block /3 accessed in the sequential execution oi a, P a will need to access 
at most 4 blocks. ■ 

Corollary 3.2. The block misalignment in the parallel execution of SPMS increases the cache miss 
bound by just a constant factor. 

Proof. By Lemma 13.11 if the sequential execution uses a cache of size M and the parallel execution 
uses a cache of size 4 M, then the parallel execution will incur at most 4 times as many cache misses. 

By regularity m, as the cache miss bounds just use the assumption that M > B 2 and are 
polynomial in B, if the parallel execution actually had a cache of size M, it would increase the 
cache miss bound by just a constant factor. ■ 


3.4.2 Caching Costs Under Work Stealing 

Suppose a parallel execution incurs S steals and let R(S) be a bound on the number of additional 
cache misses this computation incurs as a function of S. There is a simple upper bound on R(S) 
when using a work-stealing scheduler. Consider the sequence of steps a executed in a sequential 
execution. Now consider a parallel execution that incurs S steals. Partition a into contiguous 
portions so that in the parallel execution each portion is executed in its sequential order on a 
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single processor. Then, each processor can regain the state of the cache in the sequential execution 
once it has accessed 0(M/B ) distinct blocks during its execution (see, e.g., [T] in conjunction with 
Lemma [3.ip . Thus if there are K portions, then there will be R(S ) = 0(K ■M/B ) additional cache 
misses. Furthermore, it is shown in p] that I\ < 2S + 1 for a work-stealing scheduler. We review 
this bound in Lemma 13.41 using the notion of task kernels. 

The following fact is well-known (see, e.g., PI)- 

The Steal Path Fact (for Work-stealing Schedulers). Let r be either the original task or a 
stolen subtask. Suppose that r incurs steals of subtasks 7i, ■ ■ ■ , Tfc. Then there exists a path P T in 
r’s computation dag from its root to its final node such that the parent of every stolen task r* lies 
on P T , and every off-path right child of a fork node on P is the start node for a stolen subtask. 
Comment. The above fact follows by observing the nature of a depth-first-search based sequential 
execution of recursive computations, and the dequeue mode of steals under work-stealing. 

Tasks and Task Kernels. Given a parallel execution that incurs steals, the following definition 
gives a partition of the computation dag into a collection of task kernels induced by the steals. 

Definition 3.3. (Task kernels under work-stealing.) Consider a parallel execution of SPMS 
under work-stealing, and let it incur S steals, numbered as <ti, 02 , • • • ,05 in the order they occur 
relative to a sequential execution. In turn, we partition the computation dag into task kernels with 
respect to the sequence Ej = (a\, 02 ,..., af) to create the collection Ci. We let Eo be the empty 
sequence and its associated partition Co be the single task kernel containing the entire computation 
dag. For each i > 1 the partition Cj+i is obtained from Ci as follows. Let r be the task kernel in 
Ci that contains the fork node Vf at which steal < 7 j + i is performed, and let Vj be the corresponding 
join node. Then, r is partitioned into the following three task kernels. 

1 . pi, the stolen subtask created by <t,; + 1 . 

2. fi 2 , the portion of r preceding the stolen subtask in the sequential execution. (This includes 
the portion of the computation descending from the left child of vj that precedes vj.) 

3. p 3 , the portion of r descendant from Vj in the computation dag, including node Vj itself. This 
is the remainder of r. 

Then, C i+ 1 = C t - {r} U {pi, p 2 , M 3 }- 

For this parallel execution of SPMS, the collection of task kernels is Cs- 

Lemma 3.4. Consider an execution of SPMS under a work stealing scheduler, and suppose it incurs 
S steals. The resulting partition into task kernels forms at most 2S + 1 kernels. Furthermore, each 
task kernel is executed on a single processor and each kernel forms a contiguous portion of work in 
the sequential execution. 

Proof. If S = 0, the entire computation is contained in one task kernel, and 1 = 2-0 + 1. Each 
additional steal partitions an existing task kernel into three new task kernels, and hence the size of 
the collection increases by two as needed. 

To see that each task kernel is executed by a single processor, it suffices to note that the only 
edges along which the executing processor can change are the edges from a fork node into the initial 


12 


node of a stolen subtask, and the edges into the join node immediately following the end of a stolen 
subtask; this describes precisely where the partitioning into kernels occurs. 

Finally, the fact that a kernel forms a contiguous portion in the sequential execution follows 
from the dfs-based sequencing of the computation dag in the sequential execution. ■ 

Corollary 3.5. Consider a parallel execution of SPMS using a work-stealing scheduler. Suppose it 
incurs S steals. Then R(S), the additional cache misses it incurs, above the cost of the sequential 
execution, is bounded by R(S) = 0(S ■ M/B). 

3.4.3 Scheduler Cache Miss Performance 

Consider a computation with span T^, whose sequential running time is Tj, and sequential cache 
miss bound is Q. Suppose a parallel execution incurs S steals and the cost of F(S) cache misses 
due to false sharing (which will be discussed in Sections 0 E]) ■ Let T s be the total time spent across 
all processors on performing successful steals and let T u be the total time spent on unsuccessful 
steals. Let b be the cost of a cache miss and s the cost of a successful steal. Finally, let I be the 
total time spent by idling processors. Since every processor is computing, or waiting on a cache 
miss or false sharing miss, or attempting a steal, successful or otherwise, or simply idling, we can 
bound the execution time T p on p processors as the following. (The steal cost s does not appear in 
the following equation since it is subsumed into the cost for steals, but it will be used in Section 0 
when we apply this bound to randomized work stealing.) 

T p = i (ri + b ■ Q + b ■ S • ^ + +b ■ F(S) + Ts + T u + l). (1) 

Equation (0) does not incorporate the overhead due to false sharing, except as an additive term 
F(S). This term will be discussed in Sections [5] and El 

As noted earlier, randomized work stealing (RWS) is the most popular form of work-stealing, 
and for this method good bounds are known for the expected number of steals [7, lj H%] . We will 
discuss this in Section 0 after we have discussed issues relating to false sharing. 

3.5 Resource Obliviousness 

We use the notion of resource obliviousness described in m- A multicore algorithm is resource 
oblivious if the algorithm does not use any machine parameters such as p, M, and B, and yet is 
analyzed to run efficiently in a parallel execution under a range of parameters. As noted in m 
resource obliviousness differs from the notion of multicore-oblivious introduced in [9] by being 
scheduler independent. 

We obtain bounds for our resource oblivious algorithms in terms of the machine parameters p, 
M , and B, the input size n, and the number of steals S'; S in turn may depend on the scheduler 
that handles the task transfers, but this is the only dependence of our algorithm and its analysis 
on the scheduling. 

For SPMS, we have presented an algorithm with an optimal bound for the sequential cache 
complexity Q(n), and another algorithm that performs 0{n log n) operations and achieves span 
^ = O (log n log log n). In Section 0 we describe an SPMS implementation that simultaneously 
achieves both of these bounds. Thus when executed by a work-stealing scheduler Equation (0) will 
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apply. This is a resource oblivious result, except that the analysis does not consider the cost of 
false sharing. 

In Sections [5H6] we discuss the cache miss overhead due to false sharing and present the final 
version of our SPMS algorithm. Its cost overhead for false sharing is the cost of 0(S ■ B ) cache 
misses. With a tall cache we have M > B 2 , thus this false sharing cost is dominated by the cache 
miss overhead of S ■ M/B shown in Corollary 13.51 

3.6 Related Work 

As mentioned earlier, if we ignore communication costs, then the multicore is an asynchronous 
PRAM. When we incorporate communication costs, our model differs from the BSP [21], LogP el 
QSM [18], and the multi-BSP [22] in two ways: it uses cache misses instead of gap and/or latency 
parameters for the communication cost, and it uses fork-join parallelism instead of bulk-synchrony. 

The PEM model [3] uses a parallel machine similar to ours: a shared memory model, where 
each processor has a private cache, and communication cost is measured in cache misses. However, 
the PEM model again uses bulk-synchrony. Further false sharing is not addressed in results known 
for PEM, nor is the notion of resource obliviousness. 

Orthogonal to our work are results on ‘multicore-oblivious’ algorithms for a multi-level caching 
hierarchy in [9] (see also [3 ED]). These results deal with a specific scheduler, and while the 
algorithms do not use any multicore parameters, the scheduler is aware of these parameters. This 
scheduler dependence on machine parameters appears to be necessary when considering a multi¬ 
level tree of caches. Our results deal with private caches only; as a result we are able to establish 
resource oblivious bounds for SPMS that are independent of the scheduler used. 

Relation to [Tl] . The SPMS algorithm was first presented in the extended abstract el and the 
basic algorithm is the same both here and in EL However, there are some differences between the 
results presented in El and the results we present in this paper, and we list the key differences: 

1. False sharing was considered in EL but the results there applied only to the case when 
space is not deallocated (as discussed in Section [5]). Here, we additionally consider space 
deallocation on execution stacks, and the resulting false sharing costs (in Section [ 6 ]) • 

2. The SPMS implementation in El was analyzed for two specific work-stealing schedulers: 
the Priority WorkStealing (PWS) scheduler [12], which is a deterministic variant of a work- 
stealing scheduler, and RWS. The PWS scheduler requires additional properties of balance in 
forked tasks (using BP and HBP computations [12! H3]). While these results are interesting 
on their own, they are orthogonal to the enhanced scheduler-independent resource oblivious 
notion in m which we use here. Hence the results in HU relating specifically to PWS are not 
included here. Instead, here we analyze more generally in terms of work stealing schedulers. 


4 Cache-Oblivious Parallel SPMS 

We would like to simultaneously achieve the 0(-j| logjyf n) cache miss bound of Lemma 12.41 and the 
0(lognloglog?r) parallel time of Lemma 12.51 To achieve this, it suffices to implement Step 1 so 
that it runs in parallel time O(logn) on an input of size n, and incurs Oin/B ) cache misses in a 
sequential execution. This follows because the forks and joins needed to instantiate the recursive 
calls in Steps 2 and 3 can be performed with height O(logr) = 0(log?r) trees. 
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We describe our parallel cache-oblivious implementation of Step 1 and its analysis in Section f4.21 
It uses four procedures, of which the Transposing Redistribution may be of broader interest. The 
other procedures are named TR Prep (short for Transposing Redistribution Preparation), Permut¬ 
ing Write, and Small Multi Merge. We detail these four procedures in Section 14.11 

4.1 Four Procedures for Step 1 
1. Transposing Redistribution. 

Input. A vector Y which we view as consisting of the sequence Fi, F 2 > ■ ■ ■ , Fife, • • •, Y r i, ■ • • Y r k of 
subvectors. These subvectors are specified with a list of the subvector lengths | Y V} | and their 
start indices Sij in the input vector, but ordered in column major order, which is the desired 
order in the output. 

In other words, the input is (Y'L), where Y is the input vector and L is the ordered list 
(|Fl|, Sll), (|Fl|> Srl)j ' ' ' 5 (|Fl |; ®rl)> j (|Ffc|j ' ' ' j (|Ffc I > &rk) ■ 


Output. 

Bound. 

Method. 


The transposed sequence Fi, Fi, • • • , Fi • ■ ■ ■, Ffc, ■ ■ • , Ffc- 
The output is computed with 0{\Y\/B + r 2 k/B) cache misses. 

First, by means of a prefix sums computation over the values | Y t] \, we determine the output 
location for each subvector Fy • Then we copy these vectors to their output locations, using a 
nested pair of loops, the outer loop being a fork-join computation over the subvectors in their 
output order, and the inner loop being a fork-join computation over the individual elements 
in each vector. The data in the second part of the input is exactly what is needed for this 
computation: namely for each subvector, where it occurs in the input, and what is its size. 

Analysis. Since the elements are accessed in scan order within each of the rk subvectors, the cache miss 
cost is \Y\/B + rk. If r < B, assuming M > B 2 , the most recently accessed block Fj, : , for 
each i, will fit in cache, and then the cost reduces to |F|/R. While if r > B, rk < r 2 k/B, 
so we can restate the cost as 0(|F|/R + r 2 k/B). This gives the desired bound. Note that 
we would have achieved a similar cache miss bound (interchanging r and k ) if we had used 
a fork-join computation for the outer loop with the subvectors in input order. We use the 
output order because that also gives good performance in the presence of false sharing, 

Notes. (i ) The second part of the input is in column major order so as to facilitate the creation of 
the output. 

(ii) Although our algorithm for Transposing Redistribution works with any input of the type 
described above, we will call it with sorted lists partitioned using a single sequence of pivots 
for all r lists. For this class of inputs, the procedure TR Prep, given below, generates the 
input for Transposing Redistribution in the format described above, when it is given as input 
the r sorted lists, and the list of pivots in sorted order. 


2 . TR Prep. This will be used to set up the input for Transposing Redistribution. 

Input. A sequence Y of r sorted lists F, F, • ■ • Y r , and a sorted subset P = {pi,p2, ■ ■ ■ ,Pk- 1} C Y of 
pivots. It is convenient to add the dummy items po = —00 and pk = +00 to P. In addition, 
each item e in P, other than po and pk, will have a pair of straddling items in each list F- Our 
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cache miss bound will depend on the parameter d, where the straddling items are at most d 
positions apart in every Y t . In our CO-Parallel-SPMS algorithm (given in Section 14.211 . d will 
take on two different values: r in Step la(ii)I, and W 2_1 in Step la(ii)III. The sorted pivot 
sequence P, together with the straddling items, are generated in CO-Parallel-SPMS using 
the fourth procedure given below (Small Mulit-Merge). 

Output. Let Sij be the rank in Y) of the smallest item greater than or equal to Pj-i, which for brevity 
we call Pj-i’s rank in Yi. TR Prep in effect partitions Yi into the sequence of sorted sublists 
Yu , Yi 2 , ..., Yik , where the items in Y^ lie in the range [pj-\,pj). This is done by computing 
the following output: (|Yu|,sn), (|* 2 i|, S 21 ), • • •, (|^ri|, s r i), • ■ ■, (|^ifc|, sifc), - - •, (I Y rk \,s r k)- To¬ 
gether with the list Y this provides the input for Transposing Redistribution. 

Bound. Our method for TR Prep incurs 0(kdr/B + At 2 /B + y/B ) caches misses, where y = |Y"|. 

Method. (i) For each pj, 1 < j < k — 1, perform a binary search in the interval spanned by its strad¬ 
dling items in each of the r lists in order to find pf s rank in each these lists. These are the val¬ 
ues for Si,j+i, 1 < * < r, computed in the desired output order sn, S21, ■ ■ ■ , s r i, ■ ■ ■ , s 1^, • • • s r k- 

(ii) For each pj, for each Yi, compute the length of the sublist of Yi straddled by pj and Pj +\: 
this is computed by means of two lockstep scans of the results from (i), comparing the ranks 
for successive items in P. This give the r-tuple of values |y) J+ i|, 1 < i < r. 

Yij | values have been computed in the desired output order in 


(m 


Since both the and 


steps (i) and (ii), the output can now be written in a lock-step scan of these two lists. 

Analysis. Let y = |T|, and recall that |P| = k — 1. 

Step (i) incurs 0(kr log \d/B~\) = 0(kr ■ d/B + At) cache misses. If d > B this is 0(krd/B ) 
cache misses. If d < B and r > B this is 0(kr 2 /B ) cache misses. If d < B and r < B, with 
an ideal cache the most recently searched block for each list can be kept in cache, and as P is 
in sorted order, the cache miss cost reduces to y/B. Steps (ii) and (in) incur the scan bound 
of 0(r ■ k/B) cache misses. 

Thus, this is always bounded by 0(kdr/B + kr 2 / B + y/B) cache misses. 

3. Permuting Writes. 

Input. A permutation array P of size x, and an input array A array of size £x, for some integer 
£ > 1, where only elements at positions £ ■ i, 1 < i < x, are relevant. 

Output. An output array C of size £'x, with C\£! ■ P[i]] containing the element A[£i\, for 1 < i < x. In 
other words, every £-th item in A is copied to every £'-th location in array C in the permuted 
order given by P. As with the input, only elements at positions £' ■ i, 1 < i < x, in C are 
relevant. 

Bound. This is computed with 0(£x/B + £!x/B + x 2 /B) cache misses. 

Note. If £ = £! = 1, then every position of arrays A and C is relevant, and C needs to contain the 
elements in A permuted according to P. We will use £ = 1, £' > 1 in our algorithm for the 
Small Multi Merge procedure, defined below. 
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Method. We compute the output of the permuting writes in an array C' of size x with a fork-join 
computation where the leaves correspond to the elements in vector P in their input order. 
Then, with a scan, we spread the output elements out to positions £' ■ i in the output array 
C. 

Analysis. Reading the input array A incurs 0(£x/B) cache misses. The writes into array C' could incur 
up to x cache misses, while spreading the elements in the output locations in array C incurs 
0(1'x/B) cache misses. Thus the cache bound is 0(£x/B + £'x/B) plus the cache miss cost 
Z for the permuting writes into array C'. 

If array C' fits in cache, i.e. if x < B 2 , then each block in array C' is read into cache only once 
and then Z = 0(x/B). Otherwise, x > B 2 , hence y/x/B > 1, and then we obtain x < x 3 ^ 2 /B. 
Thus the cache miss cost for permuting writes is bounded by 0{lxfB + £!xjB + x 3 ^ 2 /B) = 
0(£x/B + £'x/B + x 2 /B) cache misses. 

4. Small Multi Merge. 

Input. The input P =< W±,--- , Wy, > is an ordered sequence of multi-way merging problems 
(similar to the input to Step 3 of SPMS in Section [TJ) . Each multi-way merging problem will 
contain r lists of total length at most x. and across all problems the input P will have size p. 

Output. The multi-way merging problems will be output in sequence, with the elements of each W % 
being output in sorted order; consequently, the whole output will be in sorted order. Addi¬ 
tionally, for each item e in each Wj, the ranks of the two successive items straddling e in each 
of the r input lists forming the subproblem Wi that contains e are computed. 

In other words, each of the p items in the sorted output is listed along with the ranks for its 
r pairs of straddling items, resulting in an output of length p • (2 r + 1 ). 

Bound. This is computed with 0((x + r) ■ p/B + p ■ ry/x/B) cache misses. 

Method. Let W denote a generic Wi, and let X\, ■ ■ ■ X r be the r sorted lists in W. 

(z) For each multi-way merge problem W, for each item e in W, perform r binary searches, 
one in each of the r lists Xj forming W, to find its rank in the list. 

Compute the sum of these r ranks to obtain e’s rank in W. Note that if an item e has rank 
rj in list Xj then its two straddling items in Xj are the items with ranks r ? and rj + 1. 

(ii.) Reorder the items in W to be in sorted order, keeping the results of its r searches with 
each item. To this end, using a Permuting Write with £ = 1 and £! = 2r + 1, the items are 
copied to their new locations, which are distance 2r + l apart; then, for each item, its r search 
results are copied. 

Analysis. Step (i): If the input for each merging problem fits in cache, i.e. if x < B 2 , this step incurs 
0(r ■ p/B) cache misses (as the cost of writing the output dominates). Otherwise, let Xj 
be the length of the j-th list in problem W, and let w = \W\ = Xp=i x j- Step (i) incurs 
0(wJ2j l°g \xj/B\) = 0(w 2 /B-\-w-r) cache misses. Summed over all the multi-way merging 
problems, this totals 0{x ■ p/ B + p-r) cache misses. Since x > B 2 , we have p-r < p-ryfx/B. 
So the cost is always bounded by 0(x ■ p/B + p • ry/x/B). 

In Step {ii), the cost of the permuting write for subproblem W is 0(wr/B + w 2 /B). Summed 
over all the subproblems, this totals 0 (p • r/B + px/B). The final set of writes incurs 
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0(p ■ r/B + p) cache misses (the second term arises due to the reading of each block of 
r ranks). If each subproblem fits in cache, i.e. if x ■ r < B 2 , then the cost is 0(p ■ r/B), 
and otherwise x ■ r > B 2 and then p < pJ~xr/B). So the cost is always bounded by 
0(p-r/B + P VTrlB). 

So the overall cache miss cost is 0((x + r)p/B + p ■ ry/x/B). 

4.2 Details of Step 1 in SPMS 

We now give the efficient parallel cache-oblivious algorithm, assuming the bounds stated for the 
four procedures defined in Section 14.11 We follow the structure of the algorithm in Section 12.11 but 
Step la (ii) is changed substantially. 

CO-Parallel-SPMS (Step 1) 

Step la. Form and sort the sample set S. 

(i) For each of the r sorted lists Lj, form the sublist SLj comprising every r4 c / 2 ) _1 -th item in Lj. 

The set S comprises the union of the items in the SLj lists. 

(ii) Output the set S in sorted order in a three-step process as follows. 

I Form and sort the subarray SS comprising every r-th item in the (unsorted) S using 
a Small Multi Way Merge. Note that for each e in SS each pair of straddling items 
returned by Small Multi Merge are r apart in the corresponding list forming S. 

II Partition S about SS using a TR Prep followed by a Transposing Redistribution. 

Ill Sort the resulting subsets of S, which all comprise r lists each of at most r items, using 
a Small Multi Merge. Note that for each e in S , each pair of straddling items are 1 apart 
in the corresponding list forming S, and hence W 2_1 apart in the corresponding input 
list. 

Step lb. Form the pivot set P and use P to partition A. 

(i) Form the sorted pivot set P consisting of every 2r-th item in the sorted S. 

(ii) Partition the r lists about P by means of a TR Prep followed by a Transposing Redistribution. 

Lemma 4.1. The sequential execution of the above version of Step 1 incurs 0(n/B) cache misses, 
if c > 6 and M > B 2 . 

Proof. Steps la(i) and lb(i) are simple scans that incur 0(n/B) cache misses. The scan used for 
forming the subarray SS in Step la (ii) is similarly bounded. 

In Step la(ii) the Small Multi Merge for sorting SS in step I has only one collection of r sorted 
lists, and hence has parameters x = p = n/W 2 . Therefore it incurs 0(n 2 / (r c B) + n/ (r c / 2_1 B) + 
n 3/2/ r (3c/4-l£)) = o(n/B) if c > 4. 

In step II, the TR Prep to partition S about SS has parameters d = r, k = n/r c / 2 , and 
y = n/r C//2_1 . So it incurs 0(n/(r c ! 2 ~ 2 B)) = 0(n/B) cache misses if c > 4. The Transposing 
Redistribution that follows has |A| = n/r C//2_1 and k = n/r c ^ 2 . Thus it incurs 0(n/(r c / 2 ~ 1 B) + 
n/r c / 2 ~ 2 B )) = 0(n/B) cache misses, if c > 4. 
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In step III, the second Small Multi Merge that sorts the subsets of S has parameters x = r 2 , 
p = n/r c / 2_1 . Therefore it incurs 0(n/(W 2 ~ 3 .B) + n/ (r c / 2 ~ 2 B) + n/(r c / 2 ~ 3 B)) = 0{n/B) if c > 6. 

In Step lb (ii), TR Prep has parameters d = W 2_1 , k = n/r c / 2 , and y = n. Thus it incurs 
0(n/B + n/(r c / 2 ~ 2 B)) = 0(n/B ) cache misses, if c > 4. Finally, the second Transposing Redistri¬ 
bution has |A| = n and k = n/W 2 . Thus it incurs 0(n/B + n/r c / 2 ~ 2 /B) = 0(n/B ) cache misses, 
if c > 4. ■ 

Theorem 4.1. There is an implementation of the SPMS algorithm that performs 0(n log n) op¬ 
erations, runs in parallel time 0(lognloglog?r), and incurs 0( j1 1°^ + ^ ■ S) cache misses in a 
parallel execution with S steals, assuming that M > B 2 . 

Proof. We will use the implementation of Step 1 given above. The cache miss bound for Step 1 
is 0(n/B ) by Lemma 14.11 Each substep in Step 1 is a fork-join computation which begins with a 
height O(logn) fork tree, and is followed by the complementary join computation, with size 0(1) 
or O(logn) computations at the intermediate leaves (the size O(logn) leaf computations are for 
the binary searches). Thus the span (i.e., parallel time) for Step 1 is O(logn). An 0{n) operation 
count bound should be immediate. 

The bound of O(-g ^pj ) on the cache misses now follows from Lemma 12.41 and the bounds 
of O(lognloglogn) on the parallel time and 0(nlog n) on the operation count from Lemma 12.51 
Finally, the additive term of 0(MS/B) follows from Corollary 13.51 ■ 

Comment. Theorem 14.11 assumes a linear space implementation of SPMS. (Such an implementa¬ 
tion would use two collections of size 0(n ) arrays, and at successive recursive levels would alternate 
between which arrays are used for inputs and for outputs, analogous to what is done in the linear 
space implementation of the standard merge sort.) Otherwise, if space was allocated dynamically 
as needed for each recursive procedure, the algorithm would use 0(nlogn) space, and it would 
only be for problems of size 0(M/ log M) that their whole computation would fit in cache. In 
fact, this would not affect the complexity bounds as the term 0(j| /(( 'fu ) would be replaced by 

Wi( n log n \ _ p/ n logn \ 

^yBlog(M/\ogM)i ^^BlogiV/c 

Comparison with the partitioning process in [6]. The sample sort algorithm in Blelloch 

et al. [6] (which achieves 0(log 2 n) time span deterministically, and 0(log L5 ?r) span randomized) 
also performs a partitioning following the sorting of the sample. It used a generalization of matrix 
transpose to carry out this process; it plays the same role as our transposing redistribution, however 
neither can be used in place of the other due to the different data sizes. In [6], there are y/n sorted 
subsets each of size yfn to partition about yfn items, which contrasts with the r = n 1//6 sorted 
subsets of total size n being partitioned about r 3 = yfn items. In turns out we could apply 3 
successive iterations of their procedure in lieu of the transposing redistribution, but this does not 
appear particularly advantageous. Further, the method in [6] can incur significant false sharing 
misses in the worst case, in contrast to our method, as discussed in Sections EM 

5 Bounding False Sharing Costs in SPMS 

In this section, we address the costs incurred by SPMS due to false sharing , which we will call fs 
misses , and will measure in units of cache miss cost. Recall that we are working in a multicore 
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setting with p cores, where each core has a private cache of size M, with data organized in blocks 
of B words each. When a core needs to access a data item x, this is a unit-cost local computation 
if the block containing x is in cache, but if the block is not in cache, the core needs to bring in 
the data item from main memory at the cost of a cache miss. This is the traditional caching 
communication cost that we have considered in the previous sections, as have all other papers on 
multicore algorithms (with the exception of | d3l [14] ). 

Consider the case when cores C and C' both read the same block (3 into their respective caches 
in order to access x € f3 and x' £ (3 respectively. Suppose C now changes the value in x with a 
write. Then, the local view of (3 differs in the caches of C and C 1 , and this difference needs to 
be resolved in some manner. The current practice on most machines is to use a cache coherence 
protocol, which invalidates the outdated copy of /? in the cache of C'. So, if C' needs to read or 
write any word in f3, it will incur a cache miss in order to bring in the updated copy of (3 into 
its cache. One can readily design memory access patterns by two or more cores and associated 
asynchronous delays that result in very large cache miss costs at cores due to invalidations of local 
blocks caused by the cache coherence protocol. Under such circumstances, the block structured 
storage causes larger delays than that incurred by reading in individual words at the cost of a cache 
miss each. This is considered acceptable in practice because false sharing occurs fairly infrequently. 
It is nevertheless considered to be a fairly serious problem. 

Recently, in m, systematic algorithm design techniques were developed that help reduce the 
cost of false sharing, as were techniques specific to the class of balanced parallel (BP) and hierarchical 
balanced parallel (HBP) computations. All of the computations in Step 1 of SPMS are essentially BP 
computations, and the overall SPMS is essentially a ‘Type 2’ HBP computation, so the results in [43] 
apply for the most part. However, some of the Step 1 computations (sparse permuting and parallel 
binary searches) do not exactly satisfy the definition of BP, and since the recursive subproblems in 
SPMS do not have approximately the same size (they are only guaranteed to be sufficiently small), 
the overall SPMS does not exactly satisfy the HBP definition. So we need to modify the theorems 
in m- The modifications needed are small, but as presenting the modifications would need the 
necessary background to be provided, we have chosen instead to provide self-contained proofs of the 
results for SPMS. The main result we will establish is that the overall extra communication cost 
incurred by the parallel tasks due to false sharing in an SPMS execution is bounded by the cost of 
0(5 • B) cache misses, where B is the size of a block, and 5 is the number of steals (i.e., 5 + 1 is 
the number of parallel tasks executed). In other words each parallel task incurs an amortized cost 
of at most 0(B) cache misses due to false sharing in any parallel execution of SPMS. 

The rest of this section is organized as follows. In Section 15. II we provide some basic background 
on false sharing and results from m- Then in Section f5.21 we bound the cost of fs misses in SPMS 
assuming space cannot be deallocated. Later, in Section [6] we re-establish the same bound for the 
case when space can be deallocated and then reused on execution stacks. 

5.1 Background on False Sharing 

An fs-miss occurs if a processor P had a block (3 in its cache and another processor P' writes a 
location in (3. Then in order to access f3 subsequently, P will need to reload /3, incurring an fs 
miss. If there is a sequence of k processors writing to (3 and each one in turn prevents P from 
re-accessing (3, then P will be said to incur k fs misses. It may be that such a sequence can have 
non-consecutive repetitions, e.g. Pf, p>, Pi each blocking P’s access in turn. We assume that when 
an updated (3 is accessed by another processor P, the time it takes to move (3 to P’s cache is the 
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time of b units for a cache miss (or for 0(1) cache misses). 

We will use the following definition of block delay from [13] to bound the delay due to fs misses. 

Definition 5.1. 113 '} Suppose that there are m > 1 writes to a block f3 during a time interval 
T = [t \, £ 2 ] ■ Then m is defined to be the block delay incurred by (3 during T. 

Consider a processor P that has block ft in its cache during time interval T. If no other 
processor writes to fi during T then P will incur no fs miss on (3. On the other hand, any time 
another processor writes to (3 and P accesses j3 after that write, P will incur an fs miss. Thus the 
block delay m bounds the worst-case number of times that an updated version of the block (3 may 
need to be moved into P's cache due to changes made to the contents of (3 through writes, and thus 
it bounds the cost of fs misses at P for accessing f3 in time interval T (in units of cache miss cost). 
This is a pessimistic upper bound, since it need not be the case that P will access (3 between every 
pair of successive writes to it. However, in our analysis, we will use the block delay m to bound 
the cost of fs misses at every processor that accesses (3 during the time interval T. As a result, the 
upper bounds that we obtain are quite robust, and should apply to most methods likely to be used 
to resolve false sharing, cache coherence or otherwise. In the above definition we include the notion 
of a time interval because, for example, we do not need to account for writes that occur after a join 
if all accesses of [3 by P occur before the join. 

We now present the definitions of L(r)-block sharing and limited access variables from [13]. 
and the new notion of buffered arrays (based on the gapping technique in [[133) 5 these provide the 
algorithmic tools that enable low false sharing costs. 

Definition 5.2. IT3j A task t of size r is L(r)-block sharing, if there are at most 0(L(r )) writable 
blocks that t can share with all other tasks that could be scheduled in parallel with r and that could 
access a location in the block. A computation has block sharing function L if every task in it is 
L-block sharing. 

The next definition specifies the access constraints satisfied by all the variables in the SPMS 
algorithm. In part ( i ) of this definition we reproduce the notion of a limited-access variable from 
m, but we introduce the new notion of a buffered array in part (ii). This leads to a more general 
notion of a limited access algorithm defined in part (in) than the one given in m by also including 
buffered arrays. 

Definition 5.3. 

i. ITS} The variable x is limited-access if it is accessed 0(1) times over the course of the algorithm. 

ii. An array A is buffered if the accesses to the array occur in two phases. In the first phase, 
each entry is accessed 0(1) times. The second phase is read-only, and if there are q accesses in the 
second phase, it is guaranteed that there are no accesses to the initial and to the final q entries in 
the array. 

in. An algorithm is limited access if its variables are either limited access or buffered (note that if 
an array is limited access it need not be buffered). 

Comment. The reason for allowing both (i ) and (ii) in a limited access algorithm is that both 
lead to the key property of a limited access computation established in Corollary 15.51 We will use 
buffered arrays to implement arrays on which Parallel-CO-SPMS performs several binary searches, 
since the variables stored in these arrays may not be limited access variables. We use the gapping 


21 


technique m to create buffered arrays. Let A be the array being accessed, and suppose there are 
at most q accesses to the array. Then, instead of declaring array A, declare array A which has a 
buffer of q empty locations at its start and end. It is the SPMS tasks that declare such arrays. In 
Parallel-CO-SPMS buffered arrays are needed only in the outputs of the recursive calls, and in the 
calls to TR Prep and Small Multi Merge in Step 1 that perform binary searches in parallel. The 
outputs of the recursive calls have size linear in their input. For the binary searches, the number 
of accesses is linear (or less than linear) in the size of the input to this call to SPMS. Thus, when 
using buffered arrays for these steps, we maintain the linear space bound for each call to SPMS. 

For now we make a simplifying assumption regarding space allocation, namely that there is no 
deallocation. This ensures that every memory location stores at most one variable. In Section f5. 21 
we analyze the limited access and block sharing features in SPMS and thereby bound the delay due 
to fs misses when there is no deallocation of blocks. Later, we will extend the analysis to handle 
the general case in which deallocation can occur. 

The following lemma is a generalization of an observation noted in m- Here we generalize 
that observation to include buffered arrays, in accordance with our enhanced definition of a limited 
access computation in Definition 15.31 This lemma shows that an algorithm with limited access, 
0(l)-block sharing, and executed in an environment that has no deallocation of blocks will have 
very manageable false sharing costs at each parallel task, even in the worst case. 

Lemma 5.4. Consider a parallel execution of a limited access computation with no deallocation of 
blocks, let t be any parallel task in this execution, let /3 be a block accessed by r. The worst case 
delay incurred by r due to fs misses in accessing (3 is the delay due to O(B) cache misses, where 
B is the size of a block. 

Proof. If ft contains only limited access variables, as it is accessed at most 0(B ) times, its block 
delay is O(B). If block ft contains only a (portion of a) buffered array A, in Phase 1 for the accesses 
to A, each location is accessed 0(1) times, so there are O(B) accesses, and the block delay in Phase 
1 is O(B). Phase 2 is read-only so there is no further block delay. 

Finally we consider the case when ft contains both a portion of one buffered array A overlapping 
one end (or portions of two buffered arrays A and A! overlapping both ends), and a collection C 
of other variables outside of A (and A')-, C may include entire buffered arrays. Note that there are 
0{B ) accesses to the variables in C. Thus we only need to consider Phase 2 accesses to A and A' 
since there are 0(1) accesses to each variable in A and A' in Phase 1. In this case, the buffered 
length q in A n ft is less than B, and q also bounds the number of accesses to A in Phase 2. An 
analogous bound applies to the accesses to A' nft. Hence there are amortized 0(1) accesses to each 
location in ft when ft is shared between portions of one or two buffered arrays, where the other 
variables are either limited access or belong to entire other buffered arrays. ■ 

Corollary 5.5. Consider a parallel execution of a limited access computation with 0{l)-block 
sharing and no deallocation of blocks, which incurs S steals. Then the worst case delay incurred 
due to fs misses in this parallel execution is the delay due to 0(S ■ B ) cache misses, where B is the 
size of a block. 

Proof. Let t be either the original task or a stolen task in this execution, and suppose that r incurs 
s steals. Let ri, ■ • • , Tj be r’s stolen subtasks. Since L{r) = 0(1), each Ti shares 0(1) blocks with 
the remainder of the computation. Thus r’s stolen subtasks between them share 0(s + 1) blocks, 
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counting each block according to the number of times it is shared. Suppose that there are k such 
tasks r, which incur si, S 2 , ■ ■ ■, s/~ steals, respectively. Note that k < S + 1. Now summing over all 
the tasks shows that there are 0(S + k) = O(S) shared blocks. Hence by Lemma 15.41 the worst 
case delay incurred by t due to fs misses is the delay due to 0(S ■ B) cache misses. ■ 

Discussion. Our cost model does not assign a cost for the invalidations that are performed by the 
cache coherence protocol of local copies of a block when an entry in the block is updated. If this 
delay is larger than a cache miss cost, then we would need to use this larger cost in place of a cache 
miss cost as the unit of cost for the block delay. A further refinement arises when a large number 
of parallel reads have been performed on some of the entries in a block (3 in which an entry has 
been updated. Then, the cache coherence protocol will need to invalidate (3 in the local caches of 
all of the parallel processors that read f3. If the invalidation cost depends on the number of caches 
at which a block needs to be invalidated, then this cost could become large. In our algorithm, this 
issue arises only when performing multiple binary searches on the same array (in TR Prep and 
Small Multi Merge). As already discussed, we control this potential cost by buffering the arrays 
undergoing these multiple searches. 

5.2 False Sharing Costs in SPMS: Limited Access and Block Sharing 

It is readily seen that all of the algorithms we described for Step 1 of SPMS are both limited 
access and have L(r) = 0(1) except for the Binary Search, Transposing Redistribution, Small 
Multi Merge, and Sparse Permuting Writes computations. 

We now discuss bounds on fs misses for the routines used in Step 1 that are not immediately 
seen as being limited access or having L(r ) = 0(1). 

1. Binary Searches. 

We assume here that each of the binary searches in Step 1 of SPMS is implemented as a limited 
access algorithm, e.g., by using the recursive implementation. In fact, as we will see subsequently, 
with a more precise analysis of where variables are stored, developed in Section [H we are able to 
dispense with this assumption. 

We call binary searches in parallel in our procedures for TR Prep and Small Multi Merge. We 
will use buffered arrays for the r sorted lists being searched in the binary searches in order to 
obtain limited access computations for these steps. This results in a limited access algorithm, and 
Corollary 15.51 will apply. 

2. Transposing Redistributions. 

Recall that we have formulated this as a log-depth fork-join computation according to a non¬ 
standard ordering of the leaves, namely according to the ordering of the writes. The result of this 
ordering is that this computation has L(r ) = 0(1), since only the two end blocks of the set of 
blocks accessed by a parallel task that is executing part of this computation can be shared with 
other parallel tasks. Had we accessed the data according to the ordering of the reads, we would 
not have obtained 0(1) block sharing. Since the variables are limited-access, Corollary 15.51 applies. 

3. Small Multi Merge. 

In the Small Multi Merge routine there is a step that copies over the 2r rank values for each 
element after the elements have been copied over with a permuting write. As in Transposing 
Redistribution, the fork-join tree for this copying step has its leaves in the output order. This 
causes the computation to have L(r) = 0(1). 
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4. Permuting Writes. 

In our implementation of Permuting Writes we have a step in which a sequence of x elements 
are written into an array C' of size x in an arbitrarily permuted order. (The other steps in 
this algorithm are implemented with scans that are 0(l)-block sharing and access limited access 
variables.) Consider the algorithm that writes x elements into an array where the writes could occur 
in an arbitrarily permuted order. This algorithm is limited access since each location has only one 
write. However, every block in the computation can be shared so the block sharing function is 
L(r) = r. This block sharing cost is high, and only gives a bound of 0(x ■ B) fs misses, which is 
excessive. 

To reduce the fs miss cost, we modify the Permuting Writes algorithm to split the writing into 
two substeps: 

Permuting Step in the Permuting Writes Algorithm. 

a. Initialize an all-zero array of size x 2 and write the i-th output element into location i ■ x. 

b. Compact the result in step a into an array of size x. 

Step b is readily achieved by a log-depth fork-join computation with L(r) = 0(1). We now 
bound the fs misses incurred by step a, which performs writes to locations i ■ x, 1 < i < x, where 
each of these x locations is written exactly once, but not in any particular order. 

We need to verify that the previous bound on cache misses applies to this modification of the 
permuting writes algorithm. We analyze it as follows. As before, reading A incurs 0(tx/B ) cache 
misses. Now the writes are made to an array C' , of size x 2 . If \C'\ < B 2 , then this incurs 0{x 2 /B) 
cache misses. Otherwise, it incurs 0(x) cache misses, and as x > B here, this is also 0(x 2 /B ) 
cache misses. The final writing to C incurs 0(x 2 /B + i'x/B) cache misses. This is a total of 
0(lx/B + i'x/B + x 2 /B ) cache misses. 

For the fs miss analysis, we only need to consider the case when x < B, since there is no false 
sharing when x > B. Let x < B, and let i be the integer satisfying x ■ i < B < x ■ (i -f 1). 
Then, at most i writes can occur within a single block. Each of these writes can incur a cost of at 
most i cache misses while waiting to get control of the block. Hence, the total cost of fs misses is 
0(x ■ i ) = 0(x ■ B/x ) = 0(B) by our choice of i. 

The above analysis establishes that for any invocation of a call to Permuting Writes, the total 
overhead due to fs misses across all parallel tasks that participate in this computation is bounded 
by the cost of 0(B) cache misses. We translate this into an overhead of 0(B) cache misses per 
parallel task in the computation across all invocations of Permuting Writes by assigning the full 
cost of fs misses in an invocation to a task stolen during this invocation. Thus, each parallel task is 
assigned the fs miss cost of at most one invocation of Permuting Writes (and if there are no steals 
in an invocation, then there is no fs miss during its execution). 

Finally, we need to consider the fs misses incurred due to interactions between blocks accessed 
for different calls to Permuting Writes. For this we observe that although the writes can occur in 
arbitrary order in Step 1 within the sparse array in a call to Permuting Writes, only the two end 
blocks of the array can be shared by other parallel tasks outside this call to Permuting Writes, 
hence the block-sharing function remains L(r) = 0(1). This leads to the following lemma. 

Lemma 5.6. Consider an execution of SPMS which is implemented with Permuting Writes. Sup¬ 
pose that there is no deallocation of space while the algorithm is running. Then the total cost 
for fs misses during the execution of all invocations of permuting writes is bounded by 0(S' ■ B) 
cache misses, where S' is the number of steals of subtasks of permuting writes that occur during the 
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execution. 


We have shown: 

Theorem 5.1. Consider a computing environment with tall caches in which there is no deallocation 
of space while an algorithm is running. There is an implementation of the SPMS algorithm in 
a parallel execution with S steals and using a work-stealing scheduler that performs 0(n log n) 
operations, runs in parallel time O(lognloglogn), incurs 0() | ^ • S) cache misses, and 

incurs a total cost due to fs misses delay bounded by 0(S ■ B) cache misses. 

The Final SPMS Algorithm. Our final SPMS algorithm is CO-Parallel-SPMS, described in 
Section [4] together with the changes described above for controlling false sharing costs. 

In the next section, we will analyze this final SPMS algorithm to show that it achieves the 
bounds in Theorem 15.11 even in a computing environment where there is deallocation of space on 
execution stacks. 

6 FS Miss Analysis with Space Deallocation 

When space for local variables in a computation is reallocated, the same location within a block on 
an execution stack could be reused to store different variables, even in a limited access computation. 
Thus the bound of 0(B) block delay we established in Section [5] for any block in a limited access 
computation does not necessarily hold. 

For instance, as observed in [33], consider a simple limited access computation r consisting of 
a single fork-join tree of height h, storing one data item at each leaf and at each forking node. 
Hence the total space used by this computation on E r is 0(h) (i.e, 0(1) space for each of up to h 
segments). Consider a block (3 on E r . Over the duration of this computation, /3 stores data for a 
subtree of r of depth d = 0(min{H, h'}), where h' is the height of the subtree, and hence @(2 d ) 
different variables will be stored in (3 during this computation. Thus, although limited access limits 
the number of different values of a variable that can be stored at a given memory location, due to 
the re-use of the same location in a block for different variables, the number of different variables 
stored on /3 can be even exponential in B. 

We note that execution under functional programming does not encounter the above issue since 
new space is allocated for each new call. For such environments, our analysis here, outside of 
usurpations, is not needed. However, using E r in a stack mode is quite common. In this section we 
carefully analyze and bound the block delay due to fs misses at a block on an execution stack of a 
task in a parallel execution of SPMS. Section [6.11 summarizes some basic facts that will be used in 
our analysis. Section [6.21 bounds the block delay at a single block on an execution stack of a task 
in a parallel execution of SPMS. 

6.1 Preliminaries 

As in Section [5] we will bound the cost of fs misses by bounding the block delay m (see Defini¬ 
tion 15.ip of any block. We begin by summarizing some basic information from m- 

(Fl) Space Allocation Property fl3j. We assume that whenever a processor requests space it 
is allocated in block sized units; naturally, the allocations to different processors are disjoint 
and entail no block sharing. 
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Comment. We will use the above assumption on how the runtime system allocates blocks 
in order to limit the range of time during which moves of a block on an execution stack 
across caches can contribute to its block delay for a particular collection of accesses to it. In 
particular, if a block used for an execution stack is deallocated and subsequently reallocated, 
accesses to the block in its two different incarnations do not contribute to each other’s block 
delay. 

(F2) Block Move Fact [13]. Let (3 be a block on the execution stack E T of a task r, and let T' 
be any subinterval of time during which r is executed. Suppose that the processors executing 
stolen subtasks of r during T' access block /3 a total of x times during T'. If there are u 
usurpations of r during T'. then (3 incurs a block delay of at most 2 x + u during T 1 , regardless 
of the number of accesses to [3 made during T' by the processor(s) executing r. 

Comment. The above fact is given in m ■ It simply observes that the processor C which 
currently “owns” r’s execution stack, that is the processor which can add and remove variables 
from the stack, can make many accesses to /3 but contributes toward its block delay only if /3 
has moved to a processor executing a stolen task, and hence needs to move back to C’s cache 
at the cost of adding 1 to its block delay. 

(F3) The Steal Path Fact (from Section [3]). Let r be either the original task or a stolen 
subtask. Suppose that r incurs steals of subtasks ri, ■ ■ ■ , t^. Then there exists a path P T in 
r’s computation dag from its root to its final node such that the parent of every stolen r* 
lies on P T , and every off-path right child of a fork node on P is the start node for a stolen 
subtask. 

Comment. This fact is reproduced from Section [3l 
The following lemma is a simple consequence of the above assertions. 

Lemma 6.1. Let r be a limited access computation that incurs some steals, and let f3 be a block 
allocated on t’s execution stack E r . Let P T be the steal path in r. Then, 

1. If (3 is used only for segments of nodes that are not on P T , then it incurs no fs misses. 

2. If (3 stores some data for one or more segments for nodes on P T , then its block delay is 
0(V + JJ), where V is the number of different variables used in the segments for nodes on P T 
stored in (3, and U is the number of usurpations ofr. 

Proof. For the first part we observe that since (3 does not store data for any segment on P T , it is 
not a shared block on the execution stack, and hence it incurs no fs misses. 

For the second part, we note that, between them, the variables in the segments for nodes on P T 
are accessed 0(V ) times. But these are the only variables that stolen subtasks can access. Thus, 
by Fact (F2), the overall block delay is 0(V + U ). ■ 

6.2 FS Misses at a Single Block on an Execution Stack 
6.2.1 FS Misses on the Execution Stack during Step 1 

Facts FI F3 and Lemma 16.11 provide us the basic tools for bounding the block delay of a block 
on an execution stack. We begin by presenting this bound for the fork-join tree computations in 
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Step 1 of SPMS; in the next section we present the bound for the overall SPMS algorithm. The 
following lemma is very similar to Lemma V.4 in m , but we have generalized it to state the result 
in terms of the total space used, allowing for non-constant space usage at all nodes, and relaxing 
the requirement of limited access for the computation at the leaves of the fork-join tree. 

Lemma 6.2. Let r be a fork-join computation, where both the internal nodes and the leaf nodes may 
use non-constant space, and let the total space used by r be E. If the accesses to the segments for 
the internal nodes of r are limited access then in any parallel execution of t, any block (3 allocated 
on t’s execution stack incurs a block delay of 0(min{.B, E}) cache misses during r’s execution. 

Proof. We will apply Lemma 16.11 so we need to bound V and U, considering only the variables 
used in segments stored on the steal path P T . We first consider V . If we consider a single path from 
the root to a descendant node in a fork-join tree, there is no space reuse on the execution stack, 
and the segments for the fork nodes along this path will occupy distinct but contiguous locations 
on the execution stack. Since the total space used by r is E, this also bounds the space used by all 
segments for the fork nodes on P r . 

It is possible that a segment for a leaf node l lies on f3, and since the accesses at a leaf node 
need not be limited access (if it corresponds to a non-recursive implementation of a binary search), 
there is no bound on the number of accesses to (3 made by the computation at l. But since the 
computation at a leaf of a fork tree is always executed by the processor that executes the task kernel 
to which the leaf belongs we have by Fact (F2) that any fs miss cost for the accesses to variables 
on the execution stack for a leaf node is absorbed into the cost for fs misses at /3 for stolen tasks. 
Hence V = 0(min{B, E}). 

We bound U as in m- Usurpations occur during the up-pass, and only during accesses to 
segments on the steal path P T . There is at most one usurpation per segment stored on (3, hence we 
can bound U = 0(mm{B,h}), where h is the height of the fork-join tree for r. Since each node 
stores at least one element for its segment, h < E, hence U = 0(min{.B, E}). The lemma now 
follows from Lemma 16.II ■ 

Lemma 6.3. Each Step 1 task r uses space E = 0(log |r|). 

Proof. Let r be a Step 1 task. The sequence of segments on its execution stack E r comprises 
0(log |t|) segments for the fork nodes on the path to any leaf node in r, followed by the segments 
for the leaf node (in the case of a recursively implemented binary search, the leaf node corresponds 
to a recursive computation of length 0(log |r|) and thus induces a sequence of 0(log |r|) segments; 
otherwise the leaf computation has length 0(1) and induces a single segment). Each segment has 
size 0(1) and thus the total size of the segments is 0(log |r|). ■ 

Each substep in Step 1 of SPMS comprises either a sequence of 0(1) size n fork-join tasks, each 
of depth O(logn), or a collection of such tasks which are performed in parallel, and which have 
combined size O(n). The initiation of these tasks is itself performed by a fork join-tree of depth 
O(logn). Each fork-join tree uses constant space for each segment, and the leaves use constant 
space, except for the binary searches which, use O(logn) space; hence E = O(logn). As there are 
0(1) substeps in Step 1, it follows that Step 1 has depth O(logn). 

Thus, we have the following corollary. 

Corollary 6.4. Let r be a task executed in Step 1 of SPMS, and let (3 be a block on its execution 
stack. Then the block delay incurred by (3 is bounded by 0(min{_B, log |r|}). 
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6.2.2 FS Misses on the Execution Stack for the Whole SPMS Algorithm 


The following lemma bounds the number of different variables used by a task r in the SPMS algo¬ 
rithm of Section [5] (i.e., the cache-oblivious, parallel version of SPMS, with the small enhancements 
made in Section [5] to reduce fs miss costs). We bound the number of variables rather than the space 
used by r because the same space may be allocated in turn for the segments of successive recursive 
calls. 

Lemma 6.5. Let r be a task in an execution of SPMS. Then, the sum of the sizes of all segments 
placed on E T , over the course of the complete execution of t, for nodes on any single path in t’s 
computation dag, is 0(|r|). 

Proof. If r is a Step 1 fork-join computation, or a subtask of a Step 1 fork-join computation then 
by Lemma [6T3l the total size of the segments on any path in E r is 0(log |r|) = 0(|r|), because the 
sequence of segments are simply allocated and then deallocated; there is no reuse of space in this 
case. 

If r is a recursive call in Step 2 or 3, the first segment for the nodes on a path in t’s computation 
dag is the segment for the variables r declares, which has size 0(|r|); this is followed by the segments 
for paths in r’s Step 1 subtasks; there are 0(1) such subtasks and they each induce a sequence of 
segments as described in the previous paragraph, which therefore have total size 0(log |r|). This is 
followed by the segments for t’s Step 2 subtask, which in turn is followed by the segments for its 
Step 3 subtask. These are identical in structure, so we will describe the segments for Step 2 only. 
First, there is path of length O(logr) in the fork tree for the 0(W 2+1 / 2 ) Step 2 recursive subtasks, 
and each node on this path has a segment of constant size. This is followed by the segments for a 
recursive task similar to t but having size 0(^/|r[). Thus the total size of the segments on E r is 
given by S(\t\) = |t| + 2S'(a/|t|) and S'(l) = 0(1) which has solution S(\t\) = 0(|t|). ■ 

Consider the execution stack E r for the original task or a recursive sorting task t. At any point 
in time E r will contain a sequence of segments for the nodes that have started execution but not 
yet completed. 

At the bottom is the segment oy for the variables declared by t. This is followed by a (possibly 
empty) sequence of segments for a series of recursive SPMS tasks, where for each recursive task 
pi there are a series of segments for the fork nodes on the path to /r followed by a segment for 
the variables declared by /r; finally, at the top there may be a sequence of segments for a Step 1 
subtask, if the computation is currently executing Step 1 in some recursive SPMS subtask. 

We now bound the block delay of a block f3 on E r . If t incurs no steal then /3 incurs no fs 
misses. If t incurs at least one steal then, by the Steal Path Fact (F3), there is a unique steal path 
P T such that every stolen task is the right child of a node on P T . Further, the accesses that each 
such stolen task makes to E r are to the segment corresponding to its parent node, which lies on 
P T . Thus fs misses incurred on the execution stack will all correspond to accesses to blocks that 
contain segments for parents of stolen nodes; all of these parent nodes will lie on P T . Thus, by the 
Block Delay Fact (F2), it suffices to bound the accesses to (5 for segments placed in it that lie on 
P T . We use these observations in our proof of the following lemma. 

Lemma 6.6. Let /3 be a block on the execution stack E r of a task r in a parallel execution of 
SPMS. Then the block delay of fd is bounded by 0(min{.B, |t|}). 
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Proof. By Lemma 16.11 it suffices to consider accesses to (3 while it stores some element in a segment 
on P T . If (3 is currently of this form, then the bottommost segment it stores, call it a, must he 
on P T . since P T starts at the root of the computation dag for r. Let p be the subcomputation 
of t which declared the variables stored in segment a. This instantiation of block (3 will remain 
allocated until p completes its computation. We will now bound the number of accesses to (3 on P T 
during p's execution. By the Block Move Fact (F2), we only need to consider accesses to segments 
on P T by stolen tasks, and then add in U, where U is the number of usurpations of r that occur in 
this interval. 

If p is a Step 1 fork-join computation, then by Lemma 16.21 the block delay for (3 during the 
duration of p is 0(min{.B,log |p|}), and we are done (since (3 will be deallocated once p completes 
its computation). Now suppose p is a recursive call in Step 2 or Step 3, or a collection of such 
recursive calls initiated at a fork node. Let P p = pH P T . 

If \p\ = O(B) then by Lemma 16.51 the sum of the sizes of the segments for the computations 
for nodes on P p = p n P T is also O(B), hence 0(B) accesses to f3 by stolen nodes occur before /f is 
deallocated, and the block delay for /3 is O(B). 

If |/o| = fl(H), then the following amount of data for nodes on P p is placed on f3: 0(min{log \p\,B}) 
data for the nodes forking the recursive calls (in the case that p consists of multiple recursive calls 
initiated at a fork node in some Step 2 or Step 3 computation); 0(B) data for each of the constant 
number of fork-join computations in Step 1 of p, followed by 0(min{log \p\,B}) for the forking of 
recursive subtasks in Step 2, followed by the segments for a recursive SPMS task in this same Step 
2, followed by a similar sequence of segments for Step 3. We analyze the Step 2 cost closely; the 
same analysis applies to the cost of Step 3. If the segments for the fork-tree for the recursive tasks 
fill up (3, then the segments due to the recursive task he outside of f3, and so the total contribution 
of Step 2 is an additional O(B) data. If the fork-join tree does not fill up /5, we will also have 
to handle the size 0(\/\p\) recursive SPMS call on P T . Here again, there are two cases: One is 
that the segment for the recursive subproblem fills up (3] in this case the overall amount of data is 
0(B). Finally, if the segment for the recursive call does not fill up /3, then it has size less than B, 
and hence, by Lemma 1631 the total size of segments in /3 fl P p during Step 2 is O(B). Thus, O(B) 
variables are stored in /3 and hence the V component of j3's block delay is 0(B). 

To complete the proof we need to bound U, the number of usurpations that are followed by an 
access to /3, and that occur during p's duration. When a usurpation of r occurs, the processing of 
t k had already reached the join node v J for the usurping subtask. Consequently, at this time, the 
segments on E r are all for nodes on P T . Thus, when the usurping processor C' first accesses (3 it 
must either access a segment currently on P T , or it must remove the topmost segment on P T that 
overlaps /3 (in order that it can place a new segment on (3 which will be the location for its first 
access to (3). But as already shown, only O(B) data for nodes on P T is stored on (3 ; so the first 
possibility can occur only 0(B) times. And as there is only 0(B) data among all the segments 
of P T stored on /?, there can be at most 0(B) such segments, and hence the second possibility 
can also occur only O(B) times. Thus, U = 0(B). Hence, the block delay of any block on E T is 
0(min{H,r}). ■ 

We can now establish Theorem o 

Proof, (of Theorem 11.11) . Because of the L(r) = 0(1) bound, each steal and usurpation results in 
a further 0(1) blocks being shared (and if a block is shared by k such stolen and usurped tasks, 
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it will be counted k times). Note that there is at most one usurpation per steal. Thus there are 
0(S ) sharings, each incurring an 0(B) delay due to fs misses by Lemma 16.61 for a total delay of 
0(S ■ B ). The other bounds remain unchanged from Theorem 15.II ■ 


7 Overall Performance 

In this section we bound the overall running time of the SPMS algorithm. First, using Equation JI]), 
the running time of SPMS executing on p processors, when there are S steals is: 
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Note that the term F(S ) = S ■ B for the false sharing cost is upper bounded by the cache miss 
overhead since we have a tall cache, and hence does not appear as a separate term. 

We can apply the above equation to any work-stealing scheduler for which we have bounds on 
the number of steals. This is the case for RWS, which we consider next. 


7.1 Performance of SPMS Under RWS 


Bounds on the number of steals under RWS have been obtained in a sequence of papers, first in 
Blumofe and Leiserson [7], considering only the number of parallel steps, then in Acar et al. [lj con¬ 
sidering the overhead of standard caching costs, and more recently in Cole and Ramachandran m 
considering the overhead of both caching and fs miss costs. We will now apply the results in [14] 
to obtain a w.h.p. bound on the running time of SPMS under RWS. 


Theorem 7.1. Suppose SPMS is executed on p processors using the RWS scheduler. Then, w.h.p. 
in n, if M = 0,(B 2 ), it runs in time 


0 [ - 

,P 


n log n + b 


n logn 
B log M 



M 

+ 

"B +S . 


, . , b 2 M 

log n ■ log log n H-- - -— 

s log B 


log n 


Proof. We will substitute appropriate values in Equation fl2]). The analysis in [14] of the RWS 
scheduler on HBP algorithms, considering both caching and fs miss costs, applies to the SPMS 
algorithm and shows that, w.h.p. in n, S = O(p-[T 00 +^B^ ! ^]) = 0(p-[log ndog log n+^R^^]). It 
also shows that, w.h.p. in n, the total time T s + T u spent on steals, both successful and unsuccessful, 
is 0(s ■ p • Too) = 0(s ■ p ■ logn • log logn). Finally, we note that I = 0 as any processor without 
a task is always trying to steal under RWS, and hence there is no idle time not already accounted 
for. ■ 


We can now establish Theorem 11.21 

Proof, (of Theorem ll.2p . The first term in Theorem 17.II represent the optimal performance in terms 
of operation count and cache misses when M = Q(B 2 ). Thus, for optimality, it suffices to have the 
second and third terms bounded by the first term. As, by assumption, s = 0(b ■ -|r), the second 
term reduces to 0(b ■ • log?r • log logn). Rearranging yields the desired bound. ■ 
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8 Conclusion 


We have presented and analyzed SPMS, a new deterministic sorting algorithm. The final version 
of SPMS in Section [5] simultaneously achieves optimal 0(n log n) sequential running time, optimal 
• log M n) sequential cache misses (assuming a tall cache as required for this bound), and a 
critical path length of 0(logn • log log n), which is almost optimal. In a parallel execution under a 
work stealing scheduler, this algorithm has been designed to have a low false sharing cost bounded 
by 0(S ■ B ), where S is the number of steals. Its traditional cache miss overhead in such a parallel 
execution is 0(S ■ ( M/B )) using well-known earlier results. These bounds are resource oblivious: 
the algorithm itself uses no machine parameters, and the dependence on M and B are determined 
only by the analysis of its performance. SPMS uses a new approach to sorting, by combining 
elements from both Mergesort and Samplesort. 

An important open question that remains is whether resource oblivious sorting can be per¬ 
formed in O(logn) critical path length, while maintaining optimal sequential cache complexity. 
The challenge here is the requirement that none of the parameters p, M and B can be used within 
the algorithm. This is a topic for further investigation. 
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